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In this work, we generalize the numerical approach to Gaudin models developed earlier by u^ to 
degenerate systems showing that their treatment is surprisingly convenient from a numerical point 
of view. In fact, high degeneracies not only reduce the number of relevant states in the Hilbert space 
by a non negligible fraction, they also allow to write the relevant equations in the form of sparse 
matrix equations. Moreover, we introduce a new inversion method based on a basis of barycentric 
polynomials which leads to a more stable and efficient root extraction which most importantly 
avoids the necessity of working with arbitrary precision. As an example we show the results of our 
procedure applied to the Richardson model on a square lattice. 



I. INTRODUCTION 



Gaudin-type models represent a fertile basis for an ex- 
act approach to the dynamics of quantum many-body 
systems. These models do not require fine-tuning of the 
Hamiltonian's parameters to satisfy integrability condi- 
tions. This is especially attractive from an applicational 
point of view since these parameters can represent phys- 
ical coupling of the model. A particular example would 
be the central spin Hamiltonian^ where the coupling 
constants between the central spin and the individual 
surrounding spins are related to the parameters of the 
Gaudin model in a simple way. We are talking about a 
non-degenerate model when all these couplings are differ- 
ent. Non-degenerate Gaudin-type models have received 
considerable attention recently. Apart from central spin 
physics they are relevant to a mesoscopic BCS^, Dicke 
modeP , Lipkin-Meshkov-Glick and many other physi- 
cal models^. However in many physical situations there 
are natural degeneracies coming from the spectrum of 
the model. Thus, in the BCS model the bare dispersion 
relation Cdi}^) of paired fermions on a d-dimensional fi- 
nite lattice defines the coupling constants of the Gaudin 
model according to a simple equation, e^i(k) = e^. On a 
regular lattice this equation may have several solutions, 
and thus the spectrum of Gaudin parameters is degen- 
erate. The degeneracy depends on the geometry of the 
lattice and is normally a number of the order of one. 

Recently we introduced a fruitful approach for solv- 
ing Bethe equations for the Gaudin models^ . The basic 
idea is to use an equivalent reformulation of the com- 
plex coupled Bethe ansatz equations in terms of an ordi- 
nary differential equation, a method known as B A/ODE 
correspondence. In the spirit of our previous work, here 
we generalize our approach to degenerate systems. It 
was previously shown^ that in certain non-equilibrium 
situations in the non-degenerate Richardson model the 
number of relevant eigenstates (those bringing a non- 
negligible contribution to the wave function) can be small 
drastically reducing the effective Hilbert space. However, 



in more general situations, a large number of eigenstates 
can be necessary therefore making an efficient solver de- 
sirable. 

Degeneracies can play an important role in this respect 
since the number of solutions to the Bethe equations is 
automatically reduced compared to a equally large non- 
degenerate system. Moreover we show that the matrix 
equations involved become super-sparse for high degen- 
eracies leading to further simplifications. This property 
allows us to improve the system size without loosing 
computational speed and study dynamics of extremely 
large systems compared to the usual Bethe ansatz solv- 
able models's size for non-degenerate systems. 

In the end we are looking for polynomial solutions to 
a second order differential equation (see eq. (25)) a pro- 
cedure which is in every regard equivalent to the Heine- 
Stieltjes polynomial approach discussed in^ ^ for exam- 
ple. Here the differential equation depends on a non- 
trivial set of parameters which we first extract from the 
solutions of an ensemble of quadratic algebraic equations 
according to the procedure described in sections |TI| to 
[V} The first one summarizes the results of our previ- 
ous papei^, while in the following we explicitly derive 
all the relevant equations for the degenerate case. The 
polynomial of interest can then be obtained using a new 
method based on the Lagrange polynomial basis which is 



explained in section [VI| and can be used for both degen- 
erate or non-degenerate Gaudin model. The techniques 
discussed are finally applied to the Richardson model on 
a square lattice [VTll 



II. 



FROM BETHE ANSATZ TO ORDINARY 
DIFFERENTIAL EQUATIONS 



This work deals with a set of quantum integrable mod- 
els which fall in the rational family of Gaudin models 



2 



defined by the following operator algebra 



[S'^{\i),S'^+\\,)]=ig- 



Xi — Ao 



[S^{X,),S^{Xj)]=0,K = x,y,z. 



To each distinct energy cj we can associate a degen- 
eracy dj and a (pseudo-)spin magnitude Sj. The co- 
efficients Aj are then a product of those two values 
Aj = Sjdj. In the remainder of this paper, we deal only 
(^1) with degenerate spin 1/2 systems and we therefore fix all 



For any realization of the rational Gaudin algebra a set of 
N commuting operators Ri can be defined, therefore al- 
lowing the construction of exactly solvable Hamiltonians 
of the form r]iRi. 

The reduced BCS Hamiltonian on a regular lattice 



H 



k,cr 



k,cr 



Ck,( 



k,k' 



is a natural example of a degenerate model built accord- 
ing to ([T]) and will be illustrated in section VII This pa- 
per's aim is to provide an efficient and numerically stable 
algorithm which allows one to find exact eigenstates of 
this particular type of models by exploiting its quantum 
integr ability through the algebraic Bet he ansatz. 

In this formalism, by defining a pseudovacuum |0), one 
can write exact eigenstates of the system in a remarkably 
compact form defined by the repeated action of a Gaudin 
operator on 



|{Ai...AM})cxnS+(Ai)|0). 



(2) 



Here S+(ii) = S^{u) + iSy{u) plays the role of a general- 
ized creation operator parametrized by a single complex 
variable u whose explicit expression in terms of the fun- 
damental operat ors de fining a particular system will be 
model dependenljU^nSl. In every case, the pseudo vac- 
uum |0) is defined as a lowest weight vector such that 
S~(ix) |0) = for any value of u which would be the 
Fock vacuum in the case of the reduced BCS Hamilto- 
nian. This set of general states become eigenstates of the 
Hamiltonian provided the rapidities are solutions to the 
set of non-linear algebraic equations collectively known as 
Bethe equations. Specializing to rational Gaudin models, 
the general form of Bethe equations we address is given 
by 



1 



A, — Ao 



with 



N 



^ B_ £ 



(3) 



(4) 



We will denote the ensemble of N distinct {ej} as 
"energy levels", g as the "coupling constant" and A/" 
as the "total number of levels" (for non-degenerate sys- 
tems JV = otherwise degeneracies make JV > N). 
While the physical nature of these parameters differs 
from model to model, in the Richardson case discussed 
here these parameters do indeed correspond respectively 
to single particle energy levels and coupling constant. 



Sj = 1/2 and have Aj = dj/2. An arbitrary spin would 
simply modify Aj in a way which is strictly equivalent 
to a larger degeneracy. One should note here that a de- 
generacy in the energy levels or a larger spin actually 
differ since in the degenerate case the solutions to the 
Bethe equations no longer form a complete basis of the 
Hilbert space. However the resulting subspace is orthog- 
onal to the rest of the Hilbert space and therefore even 
for non-equilibrium problems, provided the initial condi- 
tion only involves this subspace, its elements remain the 
only states required to study the unitary time evolution. 
This work only addresses the numerical issues associated 
with finding the particular subset of eigenstates corre- 
sponding to solutions of eq. ([3| applied to the particular 
pseudovacuum defined before. 

Solutions of (|3| are found by scanning from the starting 
point ^ = for which the solutions are known (A^ 
for the Richardson model) and going to finite g via suc- 
cessive steps. Rapidities A^ can be real or pairs of com- 
plex conjugates, the points in g at which two real ra- 
pidities form a complex pair (or vice versa) can lead to 
numerical instabilities when working directly with the A^ 
themselves, an algorithm which finds these points 

beforehand allowing one to control the equations in the 
critical regions has been proposed. Here however, in or- 
der to avoid these difficulties, we simply introduce the 
function 



M 



Hz) ^ E ■ 



Xh 



P'{z) 
P{z) 



(5) 



where P{z) = Y\lt^{z - Afc). The zeros of P{z) then 
correspond to the roots A^ of the Bethe Equations (BA). 

One can easily show that, whenever the set of rapidities 
Xk satisfies BA, the ODE 



N 



K'{z) + K\z)-Y, 



2F(A«) 

Z — Xn> 



= 



(6) 



is satisfied. 

If we restrict our investigation to non degenerate spin 
1/2 systems we would only need to solve the simple 
quadratic equations 



A' 

J 



A,; 



A, 



(7) 



ofEq. del) 



which are obtained by taking the limit z Cj 
setting Ai = 1/2 Wi and A^ = gA{ej). " 

Once a set {Aj}, solution to Eq. (I7| is found, the 
corresponding set of rapidities {Xk} can be obtained by 
solving only linear equations and running a standa rd ro ot 
finding algorithm. While in previous publication^^ it 
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was suggested to use elementary symmetric polynomi- 
als to build the linear system, in section [Vl] we propose 
to use a different polynomial basis which leads to much 
improved numerical stability. 



III. DEGENERATE SYSTEMS 

For the generalization to degenerate systems we obtain 
successive differential equations by taking the first n^^ 
derivatives of Eq. ([6|, to do that let us fix some notations 
first. For each pair (e^ , Aj) we then have to solve a system 
of equations which contains the first dj^ derivatives of ^ 
and the first d^j^ derivatives of Aj, this of course has to 
be done for all j. 

The general form of the n ^^ derivative of (|6| is derived 



n+l d'^Mz) 



in Appendix A Eq. (A12) for A 



,(n) 



A(e^-)^^^ 



reads 



k=0 ^ ^ 



+ (1 



n+l 



(8) 



The original form of BA equations is 

M 



E 



At — 



N 

T — 



di 



-Afc + -=0 (10) 
9 9 



where N is the number of distinct energies e^. 

We know already that for ^ = the roots will sim- 
ply correspond to the values of the single energy states 
A/c = which are occupied (the (pseudo-) spins which 
are flipped with respect to the pseudo- vacuum) . Lin- 
earizing the system for weak coupling by inserting the 
solution 



Afc = + 9^k 



(11) 



in (10), multiplying the obtained equation by g and tak- 
ing the limit ^ ^ 0, the BA equations become 



E 



2 dj 
Afe - Ai 



C = 



(12) 



where r is the number of roots occupying the single en- 
ergy Cj at g = 0. 

It is now possible to introduce the polynomial f{x) = 
nfc=i(^ ~ ^k) which has the property 



hm — - 



E 

i^k 



(13) 



with 



N , 



n ^ ^(n+l-/c) ' 

^/{n+l-k)\{l-ej)>' 



(9) 



The last term of (|8|) cancels for dj = n + l. Thus, given 



a set of pairs {(ei, (ii), (e2, (^2), • • • , (^fe, dk)}^ we define a 
closed set of quadratic algebraic equations by using, for 
each pair, the corresponding set of equations: i.e for dj = 

1 we solve equation S^^"^ which depends on A^^^ , for dj = 2 
we solve (f *^^\f^-^^) in terms of (A^-^\ A^"^^) and so on. 



IV. WEAK COUPLING LIMIT OF THE BA 
EQUATIONS 



Thus equation (12) can be written in the form of a 
differential equatiori(in the limit x ^ A^) 



djf -xf ^{Bcj^Qxf = 0. 



(14) 



The function F{x) = djf{x)-xf\x)^{Bej^C)xf{x) 
is a polynomial of the same degree r as f{x) and with the 
same roots, the two polynomial are thus proportional to 
each other and the coefficient of proportionality is r (see 
Therefore we have F{x) = rf{x) and the equation 
can be written as follows 

xr{x) - [d, + {Be, + C)x] f\x) + rf{x) = 0. (15) 



The known solutions of equation ( |15| for 5 = and 
C = 1 (which encompasses the central spin and Richard- 
son models) are the generalized Laguerre polynomials 
{x). Therefore, for weak coupling, the solutions 



-l-d. 



to the BA equations correspond to A^ = 
A/e roots of the Laguerre polynomials L7 



-l-d. 



- gAj. with 
(x). 



Once all the equations for Aj and its derivatives are 
set, the fact that they are non-linear requires the use of 
an iterative method. To do so, one always needs an initial 
guess which constitutes a good enough approximation to 
guarantee convergence to the desired solution of the BA 
equations. As previously mentioned, this will be achieved 
by slowly increasing the coupling constant g but to begin 
this process we need to solve the Bethe equations in the 
weak coupling limit for a general degenerate system. 



V. HOW TO PROCEED 

In order to solve the sets of quadratic equations defined 
by ^ we use a combination of Taylor expansion to gen- 
erate an approximative solution at g -\- 6g and Newton's 
method to refine this approximation. For both methods 
we need to solve a linear system of equations defined by 
the block- A/" x A/"— matrix 
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(Sii 



S 



IN 



(16) 



\Sni ••• Snn / 
in which each pair of indexes (i, j) defines a matrix 



OA 



(0) 



OA 



(di-l) 



(di-l) 



with matrix elements given by 



(17) 



dA) 



dS, 



in) 



dA 



= < 



n+l 



din\g''-P+^ 



(n-p) 



p > n + 1 
p = n -\- 1 
p < n 
p = n 

(18) 



and for i j 



in) 



OA- 



ip) 



-djul-r- 







p = 
p^O. 



(19) 



To have an idea of the structure of S consider a N 
system with dj = 3 Vj, we would have 



V 



* o^ 

* * 

0> 

0) 
0> 

0; 



'* 0> 

* 
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* * 0^ 

* * * 

. * * * y 

'* 0^ 

* 
0> 



'* 0' 

* 

* 
K 0^ 

* 

* 
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* * * 



which does contain a lot of zeros. It should be obvious 
that for very large systems with high degeneracies this 
sparse structure of 5* will play an important role with 
respect to computation time. 

The Taylor expansion allows us to take larger steps by 
taking into account the derivatives of A^-^^ with respect 
to g. For a step Sg the first guess will thus be given by 



A 



(»)/ 



I 

E 

k=0 
I 

E 

k=0 



k\ 

(M! 

k\ 



A 



(n,fe) 



(20) 



where we note 



d^A 



A 



{n,k) 



the function A.- differenti- 



ated n times with respect to Cj and k times with respect 
to g. 



The ^—derivatives in (20) can be found by recursively 



solving the linear system 



■rik) 



(21) 



_ p(0,/c) j^{d^-l,k) p(0,/c) j^{dN-l.k)\ 

and the elements of vector ns^^ are given by 



^{n,k) 



1^0 




(22) 



Since only the rhs of (22) depends on k we will only 



need to decompose matrix A once (using the usual ap- 
proaches for linear systems such as QR Decomposition, 
LU Decomposition or Inversion, ...) and solve the linear 
system for the new vector up to a fixed number of 
derivatives. After the first guess is obtained a typical 
Newton method can be used to finally get an accurate 
solution. 



VI. ROOT EXTRACTION 
A. General setup 

Having obtained solutions for the set of variables {A^} 
previously defined, it remains necessary to extract the 
actual rapidities from this set. In fact, Slavnov's deter- 
minant formula^^ which give compact representations 
for scalar products and matrix elements, is yet only de- 
fined in terms of {A^}. In simple terms, one needs to 
find the roots of the polynomial P{z) which ultimately 
correspond to a given solution of the Bethe equations (|3| . 

While this constitutes a standard root-finding problem, 
the position of the roots in the complex plane can lead to 
numerical difficulties. In previous paper the mono- 
mial expansion P{z) = "^^^q Pn{{\})z^ was used. In 
this case, the coefficients are simply given in terms of the 
elementary symmetric polynomials of the set {A^} which, 
in principle, can be found by solving a linear system of 
equations. For a non-degenerate model one would use M 
equations given by 
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P{e,)A{ej) = P'iej) 

M M 

A(e,) J2 ^JPM-m = J2 mef-'Pi 



M-m- 



(23) 



However as shown by Wilkinson's numerical studieJ^^, 
expressing a polynomial in terms of its monomial expan- 
sion can give rise to serious numerical problems. In fact, 
the evaluation of the polynomial at a point z far from 
becomes very sensible to the finite numerical precision of 
the coefficients at large powers. At the end, the numerical 
error becomes rapidly sufficient to affect even the struc- 
ture (real vs. complex conjugate pairs) of the roots. In 
order to circumvent that problem we use an alternative 
representation of the P{z) polynomial by decomposing it 
onto the basis of Lagrange polynomials just like for poly- 
nomial interpolation. Picking any grid of Nq = M + 1 
distinct points Zi and the corresponding values P{zi) one 
can exactly write 



Ng 



P{z)=i{z)Y^ 



WiP{Zi) 



Ng 



(24) 



2=1 



e{z) 

P'jz) 

e{z) 

P"{z) 



M+1 

E 



i=l 
M+1 



z - Zi • 



Ui 

(z - Zi)(z - Zi) 

M+1 

E 



■ (27) 



Provided Zi differs from every cj (so that F{zi)^G{zi) 
remain finite), Eq. (25) evaluated dit z = Zi leads to the 
linear (in the coefficients Uj) equation: 



Ui 



M+1 



y ^ + 2 y - 

^ (Zi — Zn)(Zi — Zh) ^ (Zi — Zn)(Zi — Zh) 




E 



[Zi Zj) [Zi Zj ) 



+ G{Zi)Ui = 

(28) 



For a grid point Zi = e^, one can use the simpler 
P'{ei) = A{ei)P{ei) and find the equation: 



M+1 



where we defined i{z) = ]^ (^ — ^i), the barycen- 



tric weights Wi 



Zj) 



and 



WiP{Zi). 



While Nq = M + 1 points is a minimal requirement to 
represent a general polynomial of order M, here only M 
points are in fact necessary since the z^ coefficient is 
trivially 1. 

From the Riccati ODE (|6| it is simple to show that the 
polynomial obeys the following second order ODE: 



P'\z) - F{z)P\z) + G{z)P{z) = 



(25) 



with 



N 



^ C Bz ^ 



• 1 ^ ~ 



G{z) 



MB 



N 

E 



(26) 



Independently of the degeneracies dj it becomes possi- 
ble to write a sufficient large system of linear equations 
whose solution fully specifies the polynomial P{z). In- 
deed the first two derivatives of the barycentric represen- 
tation (24) are easily shown to be given by 



M+1 

E 



Ui 



M+1 

■E 



= A{ei)ui. (29) 



For any given grid {zi}, one can therefore obtain the 
coefficients {ui} of the Lagrange basis representation sim- 
ply by solving a set of linear equations. With this set of 
coefficients we have a representation of the polynomial 
from which we can extract its zeros using the standard 
Laguerre's method with deffation (i.e. dividing by (x — A) 
whenever a root A is found). In this representation deffa- 
tion is a very simple task as well since it can be performed 
by shifting one grid point, say Zi, to the new found root 
A. In doing so we simply set Ui = WiP{X) = since the 
polynomial has a root at that point. At the other points 



we have Uj 



due to the change in the barycen- 



trie weight Wi. We can then reduce the degree of the 
polynomial by one by simply removing the Zi ^ X point 
and repeat the procedure using the remaining Nq — 1 
grid points. 

In this construction, we are free to choose the grid {zi} 
as we please and this fact can be exploited to ensure a 
numerically stable calculation. Indeed, this representa- 
tion of the polynomial is heavily sensitive to numerical 
precision only when evaluating the polynomial at points 
z which are far from any grid point Zi. Provided one can 
define a grid which is close enough to the actual roots 
^iig) ^ ^iig)^ the impact of finite numerical precision 
can be controlled. In any case, the remaining error on the 
rapidities extracted from this procedure can always be 
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corrected by refining tiie values using the original Bettie 
equations (|3|. One should understand that an arbitrary 
set of Zi cannot provide a numerically stable extraction of 
the rapidities. However, contrarily to any fixed polyno- 
mial basis (such as the monomial expansion for example), 
the Lagrange basis always makes it possible to choose an 
appropriate interpolation grid which leads to a stable re- 
sult. 



B. Choosing the grid 

Whenever one is interested in studying the system for 
a wide variety of coupling strengths, the root extraction 
should then be performed at a number of intermediate 
points in the coupling strength scan. One could then 
simply define the grid points (at g -\- Ag) as the roots 
found at the previous point Zi{g + A^) = Xi{g) which 
guarantees proximity of the grid to the actual solution 
and should always lead to numerically stable evaluation 
of the polynomial. 

However, when studying a single given value of the 
coupling strength it is prohibitively time-consuming to 
perform the root extraction at intermediate points and 
one should instead exploit our knowledge of the proper- 
ties of any given solution to the Bethe equations to define 
an appropriate set of points Zi{g) which mimics the po- 
sitions of roots to be found. 

At weak coupling, every rapidity is contained within 
the bandwidth of the energy levels [ei, eAr] with each roots 
corresponding to a given energy level. It is then a trivial 
task to define a grid which is close to the roots. In the 
particular cases treated here and shown in Figures ([2| 
and ([3| we choose two different sets of grid points. The 
first set, being the one which provides solutions to the 
BA at weak g. Here we use the linearized solution given 
by Equation (11) and substitute g 0.1, in this way 



we ensure proximity to the actual solution but at the 
same time avoid the risk of having to evaluate 



£iz) 



at 

At some point in the computation, the actual 
solution approaches the grid and the evaluation of 
leads to a loss of stability. Since this occurs at a different 
point in g depending on the state and the system size we 
simply use the original BA equations to determine how 
close the extracted roots are to the right solution. When 
a criterion on the precision is no longer met, we simply 
set the new grid to be equal to the rapidities calculated 
at the previous point keeping this grid for the remainder 
of the calculation. 

When interested only in performing the root extrac- 
tion at strong coupling one can define a adequate grid 
without having to find roots at intermediate values. In 
this strong coupling regime we know that a subset of 
Udiv rapidities will diverge while the remaining M — Udiv 
stay finite with real parts contained within the previously 
mentioned bandwidth (see Figures ([2) and ([3| for exam- 
ples). We therefore need a grid defined by M—Udiv points 



Zi within the bandwidth and Udiv points which diverge 
in order to follow the diverging roots. 

In^^ it was shown that for non-degenerate spin-^ sys- 
tems, provided M < N/2 (^^), the configuration of roots 
at zero coupling is sufficient to determine the number Udiv 
of roots that will diverge in the strong coupling limit via 
a simple procedure. It simply relies on defining contigu- 
ous blocks of t and | spins. Starting from the highest 
energy level, the first t block has size Pi, the second P2 
while the contiguous blocks of ^ separating them have 
sizes i^i,i^2, ... With Ni, the total number of blocks, we 
have: 

ridiv = [Pn, + c^N.-i - Min(PAr, + aAr,-i, i^ivj] (30) 
with the ai terms defined recursively as: 
ao = 0, 

ai = [Pi + ai-i - Mm{Pi + c^i-i, i^i)]. (31) 

In the degenerate case, the procedure remains exactly 
the same. We simply need to consider the degeneracies 
to be slightly lifted and consider any rapidities A = to 
flip the spin at the "bottom" of the corresponding set of 
degenerate energies as in Fig. [l]s example . 



o o 00 • • 00 

ei 62 



• • •• 



• 00 

£4 



€5 



(32) 



FIG. 1. Graphical representation of the state {A(^ = 0)} = 
{e2, £2, es, es, es, es, £4, £5} for a configuration with degenera- 
cies {di = 4,^2 = 4,^3 = 4, ^4 = 3,^5 — 1}. The black 
dots correspond to t (spins flipped from the pseudo- vacuum) . 
From the counting of blocks, this particular state would have 
1 diverging rapidity in the g ^ 00 limit. 

Remarkably, for a given Udiv^ independently of the set 
of Ci and their degeneracies, we kno^v^^^ that in the g 
00 limit, these rapidities tend to gLi where Li are the 
Udiv roots of the Laguerre polynomial 
While one could numerically evaluate the exact location 
of these roots L^, it is sufficient to simply define a set 
of points which covers the known support of these roots. 
Indeed, the real and imaginary parts of the complete set 
{Li} all fall within an easily defined bounded region of 
the complex plan^^^ : 

-A/" - 2ndiv + 2M < ^{Li) < -A/" + 2M - 2 

\^{Li)\ < 2^-ndiy{-M -2ndiy^2M). (33) 

One can then use Udiv points defined by real Zi equally 
spaced within the interval [ei + (— A/* — 2nc^^^ + 2M)^, ei + 
{—M^2M — 2)g] combined with M — Udiv points located 
within the bandwidth. This choice proves sufficient to 
guarantee a numerically stable evaluation of the polyno- 
mial P{z) in both regions where the zeros (and therefore 
rapidities) are to be found. One can then apply a re- 
fining procedure using the original A-dependent Bethe 
equations ([3|. 
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For any given Hamiltonian and any particular state, 
the aforementioned ideas should suffice to build an ade- 
quate grid at strong coupling. In this regime, the split- 
ting of rapidities into a diverging and finite block makes 
the choice of grid points central to the stability of the 
algorithm. In any case, it should also be possible to use 
the algorithm starting with any given grid and reuse the 
obtained solution to define a new grid. This would lead 
to an iterative process which, at first, evaluates the ra- 
pidities with a fairly large numerical error due to a poor 
choice of grid points. However, successive steps would see 
the error reduced by a better and better choice of grid. 



VII. REDUCED BCS ON A SQUARE LATTICE 

In order to demonstrate the capabilities of this method 
we now apply it in order to find a few eigenstates of a par- 
ticular degenerate Gaudin model. The two dimensional 
Hubbard model can be written in Fourier space as 

H = — 2t^^^{cos kx -\- cos ky)ci ^ 



M = 60 N = 19 A/" =121 

S8 ttg=l/30 .»SD iiq=l/30 



U_ 
1? 



/-^ fel.t fc2 



ki,k2,q 



If we restrict to fci = — ^2 = fc we can write (34|) in the 
form of a reduced BCS model, 

H = — 2t ^(cos kx + cos ky)cl ^Cj: ^ 



U_ 
1? 



E4,t'^-fc,/-£',i^fc',t' (35) 







which corresponds to the Richardson modeP^ with in- 
teraction parameter g = ^ and single particle energies 
= — 2t(cos kx + cos ky). 

In the examples shown below we set t = 1 and define cj 
for a set of points (kx^ky) = ^{ux^riy) with {fix^riy) = 
{0,...,L}. 

Figures ( |2|3[ ) show the behavior of the rapidities for 
the degenerate energies from (35) for the ground state 



and some excited states of two different sized systems 
{L = 10 and L = 15). In the L = 10 case there are N = 
19 distinct values of cj with 4 distinct degeneracies d = 
{1, 4, 8, 20}, in the L = 15 case there are N = 36 distinct 
values of cj with 2 distinct degeneracies d = {4, 8}. 

As one can see, the steps in g which can be taken 
while maintaining stability can be quite large compared 
to the rate of change of the rapidities with respect to g. 
Even in this degenerate case, this allows a rapid scan in 
the coupling constant opening the possibility of solving 
a large number of eigenstates in a reasonable amount of 
computation time. 

One should note that the step size could easily be in- 
creased with increasing g without affecting the stability 



FIG. 2. On the top the real and the imaginary parts of the 
rapidities of the ground state with respect to g, the black dots 
represent the points at which the roots are actually computed. 
Step size is Sg = 1/100 from to 0.05 and 6g = 1/30 from 
0.05 to 1. The second and third figure from the top show 
the real and imaginary parts of a small excitation and a large 
excitation respectively, the step size is 6g = 1/100 from to 
0.1 and Sg = 1/40 from 0.1 to 1 for the small excitation and 
Sg = 1/100 from to 0.05 and Sg = 1/40 from 0.05 to 1 for the 
large excitation. Below we show a focus around the critical 
region (where the rapidities meet /separate to become com- 
plex/real) for the ground state and the lowly excited state. 
At last we show the behavior of the highly excited state's ra- 
pidities on an extended region going from ^ = to ^ = 20, 
step size is = 1/30. 



in any way. In fact, the choice of the step size only de- 
pends on the behavior of the Aj functions and can in 
principle be fined-tuned in a particular implementation 
(se^for more information). 
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M = 128 N = 36 M = 256 

tSg=l/a3 .Kg) tSg=l,'50 




FIG. 3. Rapidities of the ground state and an excited state. 
Step size is 5g = 1/200 from to 0.05 and 5g = 1/50 from 
0.05 to 1. 



VIII. CONCLUSION 

In this work, we showed how one can exploit the 
ODE/BA correspondence to design an efficient numeri- 
cal approach for finding solutions to the Bethe equations 
defining eigenstates of Gaudin models when degeneracies 
are present. It ultimately results in finding solutions to a 
set of quadratic equations expressed in terms of new vari- 
ables. Moreover, it tuns out that in such an approach, 
degeneracies or equivalently (pseudo-) spins of magnitude 
larger than ^ give rise to a sparse matrix structure mak- 
ing them actually simpler to treat than an equally large 
non-degenerate system. Examples for the reduced BGS 
Hamiltonian on a square lattice were shown, demonstrat- 
ing the efficiency of the method. 

We also introduced a new numerical technique to ex- 
tract actual rapidities from our intermediate variables. 
Being based on the Lagrange polynomials barycentric in- 
terpolation it offers much better numerical stability than 
previous suggestions based on a monomial expansion. All 
in all the techniques discussed here offer a remarkably fast 
and efficient algorithm for systematically finding solution 
to the Bethe equations. 

The numerical stability and computation speed ob- 
tained via this method opens the possibility of study- 
ing non-equilibrium dynamics of such systems. Since 
it allows one to compute a large number of eigenstates, 
when possible, it could be combined with some trunca- 
tion scheme allowing one to reduce the effective Hilbert 
space (see^ for such an example) . In more general set- 
tings, it can also allow one to use a simple Monte Carlo 



approach which would sample a large part of the Hilbert 
space. 

ACKNOWLEDGEMENTS 

O. E. A. and V. G. are supported by Swiss National 
Science Foundation. A. F. is supported by the German 
Research Foundation (DFG). 

Appendix A: Derivation of 8^^^ 



Deriving the equation for the n derivative of 
^(.) = A'(.) + A2(.)-f:?^ = 0. 



q; = 1 



(Al) 



is mostly a technical task, which this appendix adresses 
by showing how to obtain the derivatives for the most 
general F{Xk). We are interested in 



N 



N 

E 



A, 



1 (ci - Afe) 



/(A;c), 



(A2) 



were we defined f{Xk) = ^{BX^ + C) 



Inserting ( |A2[ ) in Eq. ( |A1| we get 

M 



2Ai 



Z — Xry 



M N 



= 0. 



(A3) 



The ffist three terms of (|A3| are easily derived with 
the help of Leibniz relation 



(hg) 



(n) 



k=0 



(A4) 



We then focus here in the derivation of the last term. 
Defining A(z)^"^^ as the n*^ derivative of the function 

^(^) = X]/c=i z-Xk ^^^^ respect to we have 



M 



A{z) 



(n) 



(-irn!^ 



1 



k=l 



(^-Afe)"+i- 



(A5) 



The n derivative is then given by 

/ M N , \ 



EE 



M N 



(-ir-'EE 



{'ti (ei-Aa)(^-Aa) 



n+1 



(*)• (A6) 
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To write ( |A6[ ) in terms of A(ej) one has to take the hmit 
z Cj. The term in the sum can then be expanded in 
the form 



1 



(e* - K){ej - Ac) 
1 1 



n+l 



n ^ 



(A7) 



Inserting (A7) in (A6) we find 



(*) 

M N 



ivi ly r/ 1 1 



1 



+ 1 



E 



A(ei) - A(e,) 



^(g^.)(n-fe+l) 



where the sum over a was substituted with the functions 



A(ej) and their derivatives using Eq. (A5). 



One only needs now to get rid of the divergence given 
by the term z = j in the sum over i = 1, . . . , A^. Since 
we have Ci — Cj ^ S ^ the fraction can be 

J {Ci — Cj) ^ 

written in terms of a derivative of A. Therefore, writing 
A(e^ + 5) in its Taylor form, 



72+1 ^J. 



A(e. + 5)«E^A(e,)('=)=A(e,)- 



fe=0 



k=l 



(A9) 



one finds, for the term i ^ j of (**) 



(**)k 



k=l 



\(n-/c+l) 



1 



5^ (n-/c + l)! 
1 



^ £k — n—l 



(n+l) 



k=l 
n 

-E 

/c=l 



(n+l) 



where the second and the third term cancel themselves 
since they correspond to the same sum with reversed in- 
dices /c = 1, . . . , n ^ k = n, . . . , 1 and k ^ n — k -\- 1. 
The third term of ( A3 ) is finally given by 



N 



n 



A(eO - A(e,) 
A(ej)("-'=+i) 1 



(e,-e,)fc (n-fc + 1)! 



(All) 



Thus, using Leibniz equation (A4) and noting that 
/(Ac)^'"^ = V n > the n^^ derivative of (Ias]) in the 



limit z Cj reads 



^(5n6,A(6,)(^-^)+CA(6,)(^^) 



9 

N 



A(ei)-A(e,) 



E (n-fc + l)!j- (^^2) 
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